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ABSTRACT 


This  report  describes  a  digital  computer  program  which  uses  the 
physical  optics  approximation  to  calculate  the  scattering  properties  of 
convex  perfectly  conducting  targets  with  arbitrary  shape.  The  target 
shape  is  described  in  terms  of  the  coordinates  of  a  large  number  of 
points  on  the  surface  . 

The  program  handles  bistatic  as  well  as  backscattering  problems. 
The  input  data  specify  the  frequency,  the  incidence  angles  (0{,  §[)  and 
the  scattering  angles  (Bsi^s)*  In  the  CW  case,  the  output  data  give  the 
complex  elements  in  the  scattering  matrix. 

The  program  also  handles  the  pulse  case  where  the  incident  wave¬ 
form  has  a  finite  number  of  cycles  . 

Graphs  are  included  to  illustrate  typical  results  for  the  following 
target  shapes:  sphere,  spheroids,  ogive,  and  cone. 
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A  COMPUTER  PROGRAM  FOR  PH YSICAL- OPTICS 
SCATTERING  BY  CONVEX  CONDUCTING  TARGETS 


I.  INTRODUCTION 

Although  the  physical  optics  approximation  for  scattering  by  per¬ 
fectly  conducting  targets  has  well-known  limitations,  it  does  provide 
rapid  and  efficient  calculations  which,  in  many  cases,  are  reasonably 
accurate  when  the  target  is  large  in  comparison  with  the  wavelength. 

In  this  report  we  present  a  digital  computer  program  which  uses 
the  physical  optics  formulation  to  determine  the  scattering  properties 
of  perfectly  conducting  convex  targets  with  arbitrary  shape.  Graphs 
are  included  to  illustrate  typical  results  for  several  target  shapes. 

The  program  is  written  in  the  Fortran  IV  language. 


II.  TARGET  DESCRIPTION 

The  target  shape  is  described  in  terms  of  the  cartesian  coordi¬ 
nates  (x,  y,  z)  of  a  large  number  of  points  on  the  surface.  These  points 
are  selected  on  the  intersections  of  the  target  surface  and  the  planes 
z  =  Zi  ,  z  =  Z2  ,  etc.  The  coordinate  origin  is  located  in  the  interior  of 
the  target,  and  the  surface  is  approximated  by  triangular  facets  with 
vertices  at  the  given  points.  It  is  assumed  that  these  points  cover  the 
surface  with  a  density  such  that  each  facet  is  small  in  comparison  with 
the  wavelength.  Experience  indicates  that  reasonably  good  accuracy 
can  be  obtained  with  as  few  as  20  points  per  square  wavelength  on  the 
target  surface.  The  scattering  data  converge  to  the  physical  optics 
solution  with  approximately  80  points  per  square  wavelength.  Any 
further  increase  in  the  number  of  points  will  simply  raise  the  conputa- 
tional  time  . 

When  used  with  an  IBM  7094  computer,  the  program  will  handle 
up  to  1300  points  .  Magnetic  tape  or  disk  storage  could  be  utilized  to 
extend  this  capability. 

The  meter  is  used  as  the  unit  of  length  for  the  input  data. 

Punched  cards  are  used  for  input  data  with  the  computer  program. 
Figure  1  illustrates  suitable  input  data  for  a  prolate  spheroidal  target 
with  major  and  minor  axis  lengths  of  2  meters  and  1  meter,  respectively. 
The  first  input  data  card  gives  an  integer  N  which  specifies  the  number  of 
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planes  z  =  Zi,  z  =  z2,  .  .  .  z  =  zp^  employed  to  describe  the  target.  The 
first  line  in  Fig.  1  indicates  that  30  planes  are  used  for  the  spheroid. 
This  first  card  is  followed  by  N  cards  which  list  the  z  coordinates  ZI(I) 
of  these  planes  and  the  number  of  points  NP(I)  on  each  plane.  These 
cards  must  be  ordered  with  increasing  values  of  the  z  coordinate.  Thus, 
in  Fig.  1,  card  2  indicates  that  the  first  plane  has  z  =  -1  and  there  is 
only  one  point  on  this  plane.  Card  3  shows  that  the  second  plane  is  at 
z  =  -0.98525  and  has  12  points.  Finally,  card  31  shows  that  the  last 
plane  is  at  z  =  1  and  has  one  point. 

The  remaining  data  cards  (through  297)  list  the  x  and  y  coordinates 
of  all  the  points  on  the  target  surface.  These  coordinates  are  given  in 
pairs  (x,  y),  with  three  pairs  per  card.  For  example,  card  32  in  Fig.  1 
shows  that  the  point  on  the  first  plane  (z  =  -1)  is  at  (0,  0).  The  next  four 
cards  (33-36)  specify  12  points  on  the  second  plane,  and  card  297  locates 
the  point  on  the  last  plane  at  (0,  0).  The  last  cards  (298-299)  are  de¬ 
scribed  elsewhere  in  this  report. 

In  the  computer  program,  the  coordinates  x^-  and  y^  of  point  k  on 
the  plane  z  =  z^  are  denoted  by  F(K)  and  G(K)  or  FP(K)  and  GP(K).  The 
input  data  cards  must  give  these  points  in  an  ordered  manner,  progress¬ 
ing  in  the  clockwise  direction  (for  a  distant  observer  located  on  the 
negative  z  axis)  around  the  contour  on  the  plane  z  =  z{ .  Furthermore, 
the  first  given  point  on  the  plane  z  =  z[  should  be  the  one  nearest  the 
first  given  point  on  the  preceding  plane. 

The  planes  z  =  z{  usually  have  unequal  spacing.  The  spacing  be 
tween  these  planes,  and  the  number  of  points  on  each  contour,  should 
be  designed  to  cover  the  target  with  a  fairly  uniform  number  of  points 
per  unit  surface  area. 

The  next  section  develops  the  physical  -  optic  s  scattering  equations 
used  in  the  computer  program. 


III.  THE  PHYSICAL  OPTICS  FORMULATION 

If  the  transmitting  antenna  is  at  a  great  distance  from  the  target, 
it  will  illuminate  the  target  with  an  incident  field  which  is  essentially 
a  plane  wave.  In  the  CW  case  we  let  the  time  dependence  eJ00^  be  under¬ 
stood  and  represent  the  incident  electric  field  intensity  as  follows: 

(1)  E1  =  (0t  Eg  +  E$  )  ejkri  •£ 
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Fig.  1  •  Typical  input  data  for  the  computer  program. 
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where  {r\9  0^,cj){)^are  the  spherical  coordinates  of  the  transmitting 
antenna,  (£•[,  are  the  corresponding  unit  vectors,  * 


(2)  k  =  2tt/X  , 


X  denotes  the  wavelength,  and  Ey  and  are  complex  constants  .  An 
arbitrary  point  on  the  target  surface  is  assigned  the  coordinates  (r,0,<J>), 
the  unit  vectors  (r,  6,  $),  and  the  position  vector 


(3) 


r  = 


r  r 


Finally,  the  receiving  antenna  is  assigned  the  coordinates  (r^Qst^s) 
and  the  unit  vectors  (rs,  ^s^s)-  Thus  0{  and  <f>i  specify  the  incidence 
angles  and  0S  and  are  the  scattering  angles.  The  magnetic  field 
intensity  of  the  incident  plane  wave  is  given  by 

A 

A  A  Jk  ri*Z 

(4)  H1  =  -  <t>i  Eq)  ^ - 

where 

(5)  7)  =  n/JITT  . 


The  current  density  induced  on  the  illuminated  portion  of  the  target 
surface  is  approximated  as  follows: 


(6)  J  =  2n  x 


* 


The  unit  vectors  in  the  spherical  system  are  related  to  those  in  the 
rectangular  system  by  r  =  (xx  +  yy  +4z)/r,  §  sin<\>  +  ^  cos  cj) 
and  £  x  v  =  <j) . 
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where  n  denotes  the  outward  unit  normal  vector  on  the  surface.  The 
vector  potential  for  the  distant  scattered  field  is  given  by 

(7)  A  =  — ^ —  e"^krs 

4tf  rs 

At  a  great  distance  from  the  target,  the  scattered  field  is 

(8)  Es  =  -jWA-iiL  e-Jkrs  f  f  n  X  e^k  *s  ds  . 

Ztt  rg  J  J 

From  Eqs  .  (4)  and  (8), 

-jkrs  a  •  a  * 

(9)  Es  =  e  (0i  EJ)  -  4>i  Eg)X  S 

where 

(10)  S  =  0/M  f  f  fte-ik<'i+  *•>•£  ds  . 

The  distant  scattered  field  is  represented  by 

S  A  g  A  c 

(11)  ES  =  (0S  Eb  +  <j)s  Ej)  - - 

rs 

where  E§  and  Eq  denote  complex  constants.  From  Eqs.  (9)  and  (11) 
and  the  following  vector  identity, 

(12)  A  •  (B  X  C  )  =  (A  X  B)  •  C  , 

it  is  found  that 

(13)  Ey  =  (E^i  X  0S+  Ej  0S  X  00  •  S 

and 
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(14)  E$  =  (E^  X  9s  +  Ej,  $s  X  $i)  •  S 


It  is  convenient  to  define  the  CW  scattering  matrix  as  follows: 


From  Eqs  .  (13)  through  (15),  the  complex  elements  in  the  scattering 
matrix  are  given  by 


(16) 

Six 

-  (4>^  X  0S)  •  $ 

(17) 

Sl2 

=  ($s  x  eO  •  S 

(18) 

S21 

A  A 

=  (9i  X  9S)  •  S 

(19) 

S22 

=  (<J)S  X  ^  )  •  s 

It  is  convenient  to  define  a  "vector  area  function"  A(w)  as  follows: 

(20)  S=  (j/\)  A(w)  e^  dw  , 

where 

(21)  w  =  (r{  +  rs)  •  r 


It  may  be  seen  from  Eq.  (21)  that  w  represents  one  of  the  coordinates  of 
a  point  on  the  target  surface,  in  a  rectangular  coordinate  system  that  is 
rotated  in  space  with  respect  to  the  (x,  y,  z)  system.  The  w  axis  is  co- 
planar  with  r^  and  rs  and  bisects  the  angle  between  these  unit  vectors. 
The  unit  of  length  for  the  w  coordinate  generally  differs  from  that  for 
x, y ,  and  z  . 
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For  any  point  (x,  y,  z)  on  the  target  surface 


(22)  w  =  (sinB^  cos<j>i  +  sinOs  coscj)s)x  +  (sin0^  sincfq  +  sinBs  sincps)y 

+  (cos  +  cos  Bs)z  . 


Once  the  vector  area  function  has  been  calculated,  S  can  be  determined 
efficiently  from  Eq-  (20)  by  numerical  integration-  Equation  (lu)  would 
take  more  computation  time  since  it  involves  a  surface  integral  instead 
of  a  line  integral - 

Once  S  has  been  calculated,  the  elements  in  the  scattering  matrix 
are  determined  as  follows: 


(23)  Sn  =  Sx  Xn  +  Sy  Yn  +  Sz  Zn 

(24)  S12  -  Sx  X12  +  Sy  Y i2  +  Sz  Z12 

(25)  Szi  =  Sz  Z21 

(26)  S22  =SxX22  +SyY22+SzZ22 

where 

(27)  S  =  x  Sx  +  y  Sy  +  z  Sz 

(28)  Xn  =  -  cos^  sin 9 s 

(29)  Yu  =  -  sin<fq  sin9s 

(30)  Zn  =  -  (sin c)5!  sin  <)>  s  +  cos  ^i  cos  4>s)  cos  Qs 

(31)  X 12  -  cos  sin<p{  sin0s  -  sinO^  cos  0S  sin<t>s 
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(32) 

Yiz 

=  sin0j  cos0s  cos<{)s  -  cos  0^  cos  sin0s 

(33) 

Zl2 

=  (sincj)^  cos<j)g  -  cos  sin4>s)  cos  0^  cos  0 

(34) 

Z21 

=  coscf)^  sin<J)s  -  sincj)^  cos4>s 

(35) 

X22 

=  -  sin  Q{  cos  4>s 

(36) 

Y22 

=  -  sin  0^  sin  4>s 

(37) 

Z22 

=  -  (cosc^i  cos  cj)s  +  sincj)^  sin<)>s)  cos  0^ 

In  the  pulse  case,  the  incident  wave  is  considered  to  have  M 
complete  cycles  and  a  square  modulation  envelope.  M  is  assumed  to 
be  an  integer,  and  the  incident  field  is  expressed  by 

(38)  E1  =  Ee  +  E^)  sin  (u>t  +  k  •  r) 

.  [  u(u>t  +  k  •  r)  -  u(wt  +  kr^  •  £  -  u>  t  )  ] 


where  Ey  and  are  real  constants,  u(x)  denotes  the  unit  step  function, 
f is  the  carrier  frequency, 


(39)  w  =  2tt  f  , 

(40)  T  =  l/f  , 

and 

(41)  T  =  MT 
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The  scattered  field  is  given  by 


(42)  Es(rs,  t)  =  [  0S  E|(t  -  rs/c)  +  4>s  E|(t  -  rs/  c)J  — 

T  r  c 


where  c  denotes  the  speed  of  light  in  free  space. 

It  is  convenient  to  write  the  following  matrix  equation 
/Ee(t)\  /Fn(t)  Fl2(t)^  /  Ey 


(43) 


V 


E$(t)  J 


kF2i(t)  F22(t)// 


The  pulse  response  of  the  target  is  thus  defined  with  four  functions 
of  time  given  by 


(44)  F n(t)  =  (&  X  0S)  •  F(t) 


(45)  F«(t)  =  (*.  x  ^i)  •  F(t) 


(46)  F21(t)  =  (${  X  Js)  •  F(t) 


(47)  F 22(t)  =  ($s  X  Ui)  •  F(t) 


where 


(48)  F(t)  = 


MI- 


cos(wt  +  kw) 


[  u(u>t  +  kw)  -  u(^t  +  kw  -  co t  ) J  ds 


Equation  (48)  can  also  be  written  as  follows 


(49) 


F(t)  =  —  ^  A(w)  cosfwt  +  kw)  [  u(oot  +  kw)  -  u(^t  +  kw  -  co t* )]  dw 
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Equation  (49)  is  employed  in  the  computer  program  since  it 
permits  more  rapid  calculations  than  Eq.  (48). 

The  scattering  waveforms  Fu(t),  F^t),  F2i(t)  and  F22(t)  are 
calculated  at  a  discrete  set  of  equally  spaced  points  in  time.  The 
spacing  between  these  points  is  determined  by  an  integer  L,  included 
in  the  input  data,  as  follows 


(50)  AT  =  T/L 


The  first  and  last  points  coincide  with  the  initiation  and  the 
termination  of  the  pulse  response  . 

It  must  be  noted  that  the  above  formulation  for  the  pulse  response 
assumes  the  incident  plane  wave  has  linear  polarization*  However,  it 
appears  that  the  more  general  situation  could  be  programmed  quite 
readily . 

The  physical-optics  scattering  program  is  described  in  some 
detail  in  the  following  sections  . 


IV.  APPROXIMATING  THE  TARGET 
WITH  A  POLYHEDRON 

Figure  2  shows  the  first  91  statements  or  cards  in  the  computer 
program.  Card  26  reads  the  integer  N  which  specifies  the  number  of 
planes  used  in  describing  the  target  surface.  Card  29  reads  the  z 
coordinate  2,1(1)  for  each  of  these  planes  and  the  number  of  points  NP(I) 
on  each  plane.  The  dimension  statements  reserve  core  storage  for  200 
planes,  500  points  on  each  plane,  and  a  total  of  1300  points. 

Card  37  reads  the  x  and  y  coordinates  of  the  points  on  the  first 
plane,  and  card  48  reads  the  coordinates  of  the  remaining  points. 

In  addition  to  reading  the  target  description,  this  first  section 
of  the  program  sets  up  a  polyhedral  approximation  for  the  target.  This 
polyhedron  has  many  triangular  facets,  with  vertices  at  the  given  points. 
If  each  facet  is  small  in  comparison  with  the  wavelength,  the  scattering 
properties  of  the  polyhedron  will  be  nearly  the  same  as  for  the  real  tar¬ 
get  of  interest . 
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SExELuTE  I8J08  001 

% I B  JOB  MAP  Oo<2 

SIBFTC  PHYSOP  N0DEC<  0C3 

21  FORMAT (  IX  ,dE 15.7 )  004 

22  FORMAT ( 1 X.8F1 5.7 )  005 

23  FORMAT (1  HO)  006 

24  FORMAT ( l X. 1  OF  12.4 )  007 

25  FORMAT (8F9. 3 )  008 

26  F  ORM AT ( 7F 1 0 • 5 )  009 

27  F ORMAT (71  1 0  )  010 

28  FORMAT ( 1  Hi  )  Oil 

29  FORMAT ( 1 F 1 0.5.6  I  1 0 )  012 

31  F  ORMAT (IX. 7H  THI  =  »  F  1  0  •  4  » 7H  PHI  =,F10.4.7H  THS  =.Fi0.4,  013 

27H  PHS  = .F 1 0 .4/  )  014 

32  FORMAT ( 1 X.5HFMC  =.F10.4.8H  « AVE  = . F 1 5 • 7/ )  015 

33  FORMAT ( IX, 5HS1 1  =.2F15.7.7H  S12  =.2F15.7,7H  521  =»2F15.7/  016 

27H  S22  = . 2F 1 5 .7/ )  017 

34  FORMAT ( 1 X.6HEA1 1  =.Fl5.7.dH  EA12  =.F15.7.8H  £A2 1  =.F15.7,8H  EA2  018 

22  =.F15.7/>  019 

DIMENSION  AX ( 500 ) , AY ( 500 ) . AZ ( 500 )  020 

D 1  MENS  1  ON  X (2600 ) . Y (2600 ).Z(2600).XN(2600).YN(2600).ZN(2600)  021 

DIMENSION  F (500 ) .G ( 500 ) .FP ( 500 ) .OP (500 ) . Z I  (200  )  ,NP (200 )  022 

COMPLEX  SX.SY.SZ.Sl 1 .S12.S21 . S22 » QA . Qd . QC » QD » Q< » QL  023 

DATA  TP. FP1 , RAD/6. 2831853. 12.56637. .01 743329/  024 

1  J  =  0  025 

READ(5.27)N  026 

*RITE(6.27)N  027 

DO  60  I  =  1  »N  028 

60  READ(5.29 )Z I ( I ) .NP ( I )  029 

FN=N  030 

M  =  NP ( 1 >  031 

MM  =  M+  1  032 

ZA  =  Z  1(1)  033 

NC=(M+2)/3  034 

DO  80  1=1, NC  035 

<=3*1-2  036 

80  READ ( 5 »  26 )F(<)«G(K)»F(K+1  ).G(K+1  )«F(K  +  2),G(K+2)  037 

F ( MM ) =F ( 1 )  038 

G (MM ) =G ( 1 )  039 

1 F ( M.EQ. 1 )  MM= 1  040 

DO  1 50  1 1 =2 ,N  041 

MP=NP (II)  042 

MPPsMP+1  043 

ZC*Z  1(11)  044 

NC* ( MP+2 ) /3  045 

DO  90  1=1, NC  046 

<=3*1-2  047 

90  READ (5,26)FP(<)»  GP ( <  )  » FP ( <♦ 1  )«  GP (<+l  ).FP(<+2). GP ( <+2 )  048 

FP ( MPP  )  =FP (  1  )  049 

GP(MPP  )  =GP(  1  )  050 

lF(MP.EQ.l)  MPP= 1  051 

L  =  1  052 

LP=1  053 

95  X A  =  F (L  )  054 

XC=FP(LP)  055 

YA  =  G (L  )  056 

Y C  =GP ( LP  )  057 

J  =  J+1  058 

IF(L.EO.MM)  GO  TO  98  059 

1 F (LP.EQ.MPP )  GO  TO  100  060 

RA= (FP (LP* 1  >-F(L)  >**2*(GP(LP+1  >-6(L  )  )**2  061 

Fig.  2.  First  section  of  the  computer  program. 


11 


RB= (FP (LP  >-F ( L+ 1  ))**2+(GP(LP)-G(L  +  l  )  1**2  062 

IF(RA.GT.RB)  GO  TO  I  00  063 

9B  LP=LP+I  064 

XB=FP(LP)  065 

YB=GP(LP)  066 

ZB  =  ZC  067 

GO  TO  110  066 

100  L  =L ♦ 1  069 

XB=F(L)  070 

YB=G(L)  071 

ZB=ZA  072 

110  X  (  J  )  =  (  XA-f  XB  +  XC  ) /3.  073 

Y ( J ) = ( YA  +  YB+YC  )/3»  074 

Z( J)= (ZA+ZB+ZC )/3.  075 

XN<  J  )  =  < ( YB-YA  )* (ZC-Z A  )- ( Zri-Z A  )* ( YC-YA )  >*.5  076 

YN  (  J  )  =  (  (  ZB-ZA  )*  (  XC-XA  )  -  (X6-XA  )*  (ZC-ZA  )  )  *  •  5  077 

ZN ( J  )  = (  (XB-XA  )*( YC-YA  )- (YB-YA  )* (XC-XA ))*. 5  076 

IF (L.LT.MM. OR.LP.LT. MPP )  GO  TO  95  079 

DO  120  1=1, MPP  OBO 

F ( I ) =FP ( 1 )  081 

120  G ( I ) =GP ( 1 )  062 

M  =  MP  063 

MM  =  MPP  084 

ZA=ZC  085 

150  CONTINUE  066 

NT  =  J  087 

DZ= • 0  088 

DO  155  I =2  « N  089 

DEL  =  Z 1  (  I  >  — Z I  (  1-1  )  090 

155  IF (DEL.GT.DZ )DZ=DEL  091 


Fig-  2.  First  section  of  the  computer  program,  (cont.) 
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Since  each  facet  is  assumed  to  be  small,  the  only  parameters 
that  must  be  determined  and  stored  for  each  one  are  its  vector  area 
and  the  coordinates  of  its  midpoint. 

Each  triangle  is  assigned  an  index  number  denoted  by  J.  The 
coordinates  of  a  point  at  the  center  of  triangle  J  are  denoted  by  X(J), 
Y(J)  and  Z(J).  The  cartesian  components  of  the  vector  area  (using 
the  outward  normal  direction)  of  triangle  J  are  XN(J),  YN(J)  and  ZN(J). 
The  dimension  statements  reserve  core  storage  for  these  parameters 
for  2600  triangular  facets.  The  number  of  triangles  is  approximately 
twice  the  number  of  points  on  the  target  surface,  and  is  denoted  by  NT* 

The  technique  for  fitting  a  polyhedron  to  a  given  array  of  points 
is  illustrated  in  Fig.  3.  The  computer  first  processes  all  the  tri¬ 
angular  facets  on  the  first  "zone"  of  the  surface  (i.e.,  the  portion  of 


Fig.  3.  A  typical  triangular  facet  on  the  polyhedron. 
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the  surface  bounded  by  the  first  two  planes  z  =  and  z  =  z2)  •  It  then 
proceeds  to  the  next  zone.  In  the  situation  illustrated  in  Fig.  3,  the 
computer  is  identifying  and  processing  the  triangles  in  Zone  II-  1  .  The 
x  and  y  coordinates  of  the  points  on  contour  II- 1  are  denoted  by  F  and 
G,  and  those  on  contour  II  by  FP  and  GP.  The  first  triangle  on  this 
zone  has  one  vertex  at  the  point  x  =  F(l),  y  =  G(l)  and  z  =  ZI(II-l). 
Another  vertex  is  at  x  =  FP(1),  y  =  GP(1)  and  z  =  ZI(II) .  The  third 
vertex  of  this  first  triangle  is  at  (F(2)rG(2))  or  (FP(2),  GP(2))f  which¬ 
ever  yields  the  "most  compact"  triangle. 

When  the  computer  has  completed  the  processing  of  a  given 
triangle  in  a  given  zone,  two  vertices  on  the  next  triangle  are  prede¬ 
termined  since  they  coincide  with  two  of  the  vertices  on  the  last 
completed  triangle.  In  Fig.  3,  these  two  vertices  are  indicated  by 
points  A  and  C.  Point  A  is  the  L-th  point  on  contour  II- 1  and  its 
coordinates  are  XA=F(L),  YA=  G(L)  and  ZA  =  ZI(II-  1 ) .  Point  C  is  the 
LP-th  point  on  contour  II  and  its  coordinates  are  XC  =  FP(LP), 

YC  =  GP(LP)  and  ZC  =  ZI(II).  To  establish  the  third  point  on  the 
triangle,  the  computer  compares  the  lengths  of  the  lines  AD  and  BC 
and  selects  the  shorter  one  to  form  one  side  of  the  triangle.  As  shown 
in  Fig.  3,  B  denotes  the  (L+l)-th  point  on  contour  II-  1,  and  D  is  the 
(LP+l)-th  point  on  contour  II. 

For  the  situation  shown  in  Fig.  3,  BC  is  shorter  than  AD.  There¬ 
fore,  B  is  selected  as  the  third  vertex  of  the  triangle  and  assigned  the 
following  coordinates:  XB  =  F(L+1),  YB  =  G(L+1)  and  ZB  =  ZA.  The 
coordinates  of  the  midpoint  of  this  triangle  are  calculated  in  cards  73, 

74  and  75.  The  vector  area  of  the  triangle  is  given  by  (P  *  Q)/2,  where 
the  vector  P  extends  from  A  to  B  and  the  vector  Q  extends  from  A  to  C. 
This  calculation  is  arranged  by  cards  76,  77  and  7  8  to  determine  the 
rectangular  components  of  the  vector  area. 

Having  completed  the  processing  of  triangle  J,  the  computer 
proceeds  with  the  next  triangle  in  a  similar  fashion. 

When  the  computer  finishes  the  first  section  of  the  program,  it 
has  in  storage  a  complete  list  of  the  midpoint  coordinates  and  the  vector 
areas  of  all  of  the  triangular  facets  on  the  polyhedron. 


V.  COMPUTING  THE  VECTOR 
AREA  FUNCTION 

The  second  section  of  the  computer  program  reads  the  coordinates 
(0i>  4>i)  the  transmitting  antenna  and  (Qs’^s)  of  the  receiving  antenna 
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and  calculates  the  corresponding  vector  area  function  A(w)  for  the 
target.  This  part  of  the  program  is  shown  in  Fig.  4. 

Card  92  reads  0^  4>{,  0S  and<f>s*  These  angles  are  denoted  by  THI, 
PHI,  THS  and  PHS .  They  must  be  given  in  degrees  (rather  than  radians) 
in  the  input  data.  For  example,  card  298  in  Fig.  1  assigns  a  value  of  20 
degrees  for  0^  and  0S,  and  zero  degrees  for  <\>[  and  4>s  . 

In  the  DO  LOOP  in  cards  115  through  137,  the  computer  scans 
the  w  coordinates  for  all  the  illuminated  facets  on  the  polyhedron  to 
determine  the  maximum  and  minimum  values  of  w  on  the  illuminated 
portion  of  the  target  surface.  These  are  designated  WMAX  and  WMIN. 
WX,  WY  and  WZ  represent  the  coefficients  of  x,  y  and  z  in  Eq.  (22). 

Cards  138  through  147  calculate  the  coefficients  Xn,  Yu,  etc., 
defined  by  Eqs  .  (28)  through  (37). 

The  rectangular  components  of  the  vector  area  function  A(w)  are 
denoted  by  AX(K),  AY(K)  and  AZ(K).  The  program  is  designed  to 
calculate  several  values  of  AX,  AY  and  AZ,  representing  samples  of 
the  functions  for  uniformly  spaced  points  on  the  w  axis.  The  dimension 
statements  allow  for  a  maximum  of  500  sampling  points.  The  area 
function  is  nonzero  only  for  values  of  w  between  WMIN  and  WMAX. 
Therefore,  this  portion  of  the  w  axis  is  divided  into  KX  segments  and 
the  sampling  points  are  selected  at  the  centers  of  these  segments.  Thus, 
AX(1),  AY(1)  and  AZ(1)  represent  the  components  of  A(w)  at  the  center  of 
the  first  segment  (i.e.,  the  segment  which  begins  at  WMIN).  Likewise, 
AX(KX),  AY(KX)  and  AZ(KX)  represent  the  components  of  A(w)  at  the 
center  of  the  last  segment  (which  terminates  at  WMAX).  The  length  of 
each  segment  (and  the  spacing  between  the  sampling  points)  is  denoted 
by  DEL.  WI  and  WF  denote  the  values  of  w  at  the  first  and  last  sampling 
points  . 

If  the  number  of  sampling  points  is  too  small,  the  area  function 
will  not  be  represented  adequately  and  one  cannot  expect  accurate  scat¬ 
tering  data.  Therefore,  the  program  always  calculates  the  largest 
number  of  samples  consistent  with  the  input  data.  To  obtain  a  more 
detailed  area  function  for  a  given  target,  one  simply  supplies  a  more 
detailed  description  of  the  target  via  the  input  data. 

The  following  expression  for  the  vector  area  function  can  be  ob¬ 
tained  by  comparing  Eqs.  (10)  and  (20): 
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(51) 


A(w)  = 


d 

dw 


j  j  ftds  • 


The  significance  of  Eq  •  (51)  can  be  clarified  by  considering  the  equivalent 
finite-difference  approximation: 

(52)  A(K)  =  _ 1_  \  \  nds  , 

—  DEL  J  JSk 

where  A(K)  represents  the  area  function  at  the  point  w  =  and  Sk 
denotes  a  narrow  zone  on  the  illuminated  portion  of  the  target  surface 
bounded  by  the  planes  w  =  wk  -  DEL/  Z  and  w  =  wk  +  DEL/  Z  . 

For  the  polyhedron  with  small  facets,  Eq.  (5Z)  takes  the  following 

form: 


(53) 


A(K) 


x  XN(I) 


+  y  YN(I)  +  z  ZN(I) 


DEL 


The  summation  indicated  in  Eq.  (53)  is  programmed  in  cards  16Z  through 
199  in  Fig.  4.  In  Eq .  (53)  it  is  understood  that  the  summation  extends 
only  over  the  illuminated  facets  with  w  coordinates  in  the  range  wk  ± 
DEL/2. 

A  simple  shadowing  test  is  programmed  in  cards  163  and  164. 
RDTN  represents  the  quantity  *  n.  RDTN  is  calculated  for  each  facet 
on  the  polyhedron.  The  facet  is  considered  illuminated  or  shadowed 
according  as  RDTN  is  positive  or  negative.  This  test  is,  of  course, 
adequate  only  for  convex  targets. 

Figure  5  shows  the  vector  area  components  AX(w)  and  AZ(w) 
for  a  prolate  spheroid.  AY(w)  is  zero  for  this  problem.  The  input 
data  in  Fig.  1  were  used  for  these  calculations. 

Although  independent  data  are  not  available  for  comparison  with 
the  curves  in  Fig.  5,  there  are  two  indications  that  these  curves  are 
reliable.  First,  these  area  functions  lead  to  accurate  physical-optics 
scattering  data.  Second,  the  area  functions  calculated  in  the  same 
manner  show  close  agreement  with  known  results  for  the  following 
cases:  backscatter  and  bistatic  scattering  from  spheres  and  axial 
backscatter  from  spheroids,  ogives  and  cones  • 
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2  READ < 5. 26 ) TH l ,PH I . THS.PHS  092 

iKRITE(6«31  )THI  iPHI  t  THS«PHS  093 

STHI=SIN«RAO*THl )  094 

C TH I =COS ( RAD* TH I )  095 

STHS=SIN<RAD*THS)  096 

CTHS=COS <RAD*THS )  097 

SPH I =S I N <RAD*PH I )  090 

C PH  I = COS ( RAD*Ph I  )  099 

SPHS=S I N <RAD*PHS )  100 

CPHS=COS<RAD*PhS)  101 

SSCS=STHS*CPHS  102 

SSSS=STHS*SPHS  103 

S  I  C I =  STH I *CPH 1  104 

S I S I =  STH  I *SPH I  105 

WX=S1CI+SSCS  106 

WY=S1S1+SSSS  107 

WZ=CTHI+CTHS  100 

DEL=0Z*S0RT < WX**2+WY**2+WZ**2 )  109 

WMAXs-i 00^000.  110 

WM1N=1 OUOCoO.  1 1 1 

<  =  2  112 

ZA  =  Z I<1)  113 

ZC  =  Z 1  <  2  )  114 

00  200  I = 1 .NT  115 

RDTN*XN(1  )*SICI+YN<  I  )*S1SI*ZN(I  )*CTHI  116 

1 F < ROT N • LE • 0 • )  GO  TO  200  117 

ZNI =  ZN  <  1  )  110 

ZB=Z ( I )  119 

IF  <ZC»GE»ZU)GO  TO  170  120 

165  K=K+1  121 

Z A  =  Z I  <  K  —  1  )  122 

ZC  =  Z I  (K  )  123 

IF <ZC.LT. Zd )GO  TO  165  124 

170  OEN  =  X<  1  )*XN<  1  )+Y(  1  )*YN(  1  )  125 

D=DEN+ZB*ZNl  126 

F= <0-ZA*ZNl l/OEN  127 

WA=X<I)*F*WX+Y<I)*F*WY+ZA*WZ  128 

F=(0-ZC*ZN1 )/0EN  129 

*1  C  =  X  <  1  >*F*WX+Y  (  1  )*F*WY+ZC*WZ  130 

IF (WC.GE.WA)GO  TO  175  131 

A  132 

*A=WC  133 

«C>WW  134 

175  IFOtfC.GT.WMAX )WMAX=WC  135 

I F ( WA .LT 1 N )WM IN=WA  136 

200  CONTINUE  137 

X  1  1 =-CPH 1*STHS  130 

Y  1  1 =— SPH I *STHS  139 

Z 1  1 =- <  SPH I *SPHS+CPH 1 *CPHS ) *CThS  140 

X 1 2=CTH I *SPH I*STHS-STHI *CTHS*SPHS  141 

Y 1 2= STH I *CTHS*CPHS-CTH I *CPH I *STHS  142 

Z 1 2  =  < -CPH I *SPHS  +  SPH 1 *CPHS )*CTH1*CTMS  143 

Z2 1 =-SPH I *CPHS+CPH 1 *SPHS  144 

X22=-STH I tCPHS  145 

Y22  =  — STH I *SPH  S  146 

Z 22  =  —  < CPH I *CPHS  +  SPH 1 *SPmS ) *C  TH I  147 

DO  205  K= 1 .500  140 

AX ( K ) =0.0  149 

AY<K)=0.0  150 

205  AZ(<)=0.0  151 

KX= ( WMAX-WV I N )/DEL  152 

Fig.  4.  Second  section  of  the  computer  program. 
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FK  =  KX  1 33 

DEL=  (  WMAX-miM  I  N  ) /FK  1 54 

mIRI  TE  (6.22  )DZ  .DEL.  mIM  I  N.  miMAX  1  35 

DEL  T  = DEL/2 •  156 

Ml  I  =WM  I  N  +  DELT  157 

WF=WMAX-DELT  158 

<  =2  159 

ZA  =  Z  Ml)  160 

ZC=Z 1(2)  161 

DO  230  I =1 .NT  162 

RDTN=XN(I)*SICI+YN(I)*SI3I+ZN(I)*CTHI  163 

IF (RDTN.LE.O. )  GO  TO  230  164 

ZN I =  ZN( I )  165 

ZB=Z(  I  )  166 

I F (ZC.GE.Zd )G0  TO  213  167 

21 0  K  =  K+1  168 

ZA  =  Z  I  (K-l  )  169 

ZC  =  ZI  (<  )  170 

IF(ZC.LT.ZB)GO  TO  210  171 

215  DEN  =  X (  I  )*XN(  I  )+Y (  I  >*YN(  I  )  172 

D=DEN+Z8*ZN I  173 

F= ( D-ZA*ZN I )/DEN  174 

mIA  =  X(1)*F*mIX  +  Y(I)*F*miY  +  ZA*wZ  175 

F= (D-ZC*ZNl )/DEN  176 

WC  =  X(1  )*F*»X  +  Y(I  )*F*WY  +  ZC*»/Z  177 

IF(MlC.GE.WA)GO  TO  218  178 

WW=WA  179 

mIA  =  WC  180 

WC  =  MIW  181 

218  KA= ( WA-W I  )/DEL+ 1 »5  182 

<C=  (  mIC-mI  I  )/DEL+ 1  .5  183 

I F ( KA.NE.KC )G0  TO  220  184 

AX (KA ) =AX (KA ) +XN( I )  185 

AY ( KA ) =A Y <KA ) +YN ( 1 )  186 

AZ ( KA ) =AZ (KA ) +ZN< I )  187 

GO  TO  230  188 

220  FKA=KA-l  189 

W=FKA*DEL+W I +DELT  190 

D= (W-WA )/ (WC-WA )  191 

AX (KA ) =AX (KA ) +XN( I ) *D  192 

AY ( KA ) =AY (KA ) +YN ( I ) *D  193 

AZ ( KA ) =AZ (<A ) +ZN ( I ) *D  194 

F  =  1  . — D  195 

AX ( KC ) =AX (KC ) +XN ( I ) *F  196 

AY ( KC ) =AY (KC ) + YN ( I ) *F  197 

AZ(KC)=AZ(<C)+ZN( I >*F  198 

230  CONTINUE  199 

W=WI  200 

DO  235  K= 1 .KX  201 

AX ( K ) =AX (K )/DEL  202 

AY(K)=AY(K)  /  DEL  <i03 

AZ(K)=AZ(K)/DEL  204 

FK=K  205 

*RITE(6.22>  FK  ,  mi  .  AX  (K  )  .  AY  (  K  )  •  AZ  (K  )  206 

235  mi  =  m/  +  DEL  207 

AX ( KX  +  1  )  =2. *AX ( KX ) -AX (KX- I  )  208 

AY(KX+1  )  =2 • * AY ( KX ) — AY ( KX— 1  )  209 

AZ (KX  +  1  )=2»*AZ (KX >-AZ (KX-1  )  210 

m/RITE  (6.23)  211 

Fig.  4.  Second  section  of  the  computer  program,  (cont.) 
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Fig.  5.  Calculated  area  functions  for  a  prolate  spheroid 
using  input  data  from  Fig.  1  . 
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VI.  COMPUTING  THE  CW  SCATTERING  MATRIX 


The  third  section  of  the  computer  program,  shown  in  Fig.  6,  reads 
the  frequency  in  megahertz  and  calculates  the  CW  scattering  matrix  and 
the  corresponding  echo  areas.  In  addition  to  the  frequency,  card  212 
also  reads  the  integers  L  and  M  which  are  defined  in  the  next  section. 

L  and  M  are  not  needed  for  the  CW  calculations  and  could  just  as  well 
be  read  later  at  card  249.  Suitable  input  data  for  card  212  are  shown 
in  card  299  in  Fig.  1  . 

This  section  of  the  program  calculates  the  vector  S  defined  by 

iwmax 

A(w)  e->kw  dw  . 
vmin 

The  integral  in  Eq .  (54)  is  evaluated  with  the  "piecewise  linear  method" 
described  in  Reference  1  .  We  integrate  over  segment  k  on  the  w  axis 
and  then  sum  on  k.  Segment  k  has  length  DEL  and  is  centered  at  w  =  wj^. 
On  segment  k  a  straight-line  approximation  is  used  for  the  area  function: 


(55)  A(w)  =  ak  +  bk  w 


By  fitting  Eq .  (55)  to  the  stored  samples  A(K)  at  wk  and  A(K+1)  at  wk+j  > 
we  find  that 


(56)  afc  =  [  A(K)  wk+i  -  AK+1 )  wjJ  / DEL 


and 


(57)  bk  =  [A(K+1)  -  A(K)]/DEL 


When  k  is  used  as  a  subscript  it  represents  the  index  of  a  segment  on 
the  w  axis  and  is  denoted  in  the  program  by  K.  When  used  otherwise,  as 
in  Eqs  .  (54)  and  (59),  k  is  defined  by  Eq .  (2)  and  is  denoted  in  the  program 
by  TPL. 
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The  integration  in  Eq.  (54)  can  be  performed  analytically  when  A(w)  is 
given  by  Eq .  (55).  This  leads  to  the  following  expression 

KX 

(58)  S  =  _i -  >  [QKA(K)  +  QLA(K+1)]  , 

~  ZttDG  Lj  ~  ~ 

k=l 

where 


(59)  DG  =  k  DEL  , 

(60)  QK  =  (GL  -  j)  QD  -  QC  , 

(61)  QL  =  QC  -  (GK  -  j)  QD  , 

(62)  QC  =  GB  ej  GB  -  GA  e-*  GA 

_  j  GB  j  GA 

(63)  QD  =  eJ  -  eJ  , 

(64)  GA  =  k(wk  -  DEL/ 2)  , 

(65)  GB  =  k(wk  +  DEL/2)  , 

(66)  GK  =  k  wk  , 

and 

(67)  GL  =  k  wj^+i 


The  summation  in  Eq.  (58)  is  programmed  in  cards  223  through 
236.  The  rectangular  components  of  S  are  denoted  by  SX,  SY  and  SZ. 
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Equations  (23)  -  (26)  are  programmed  in  cards  238  -  241  to  cal¬ 
culate  the  complex  elements  in  the  scattering  matrix.  This  scattering 
matrix  is  the  principal  output  of  this  section  of  the  program.  With  it 
available,  one  can  use  Eq.  (15)  to  calculate  the  scattering  properties 
of  the  target  for  any  given  incident  polarization. 

The  remaining  output  from  this  section  does  not  necessarily 
have  any  significance  except  in  special  cases.  In  backscatter  problems, 
however,  EAl  1  represents  the  echo  area  in  square  meters.  Further¬ 
more,  EA1 1  and  EA22  represent  the  E-plane  and  H-plane  echo  areas 
when  =  <j>s  and  the  target  has  axial  symmetry  with  respect  to  the  z 
axis  . 


Using  the  techniques  described  above  and  the  input  data  shown 
in  Fig.  1,  the  computer  yields  the  following  CW  scattering  data  for 
the  prolate  spheroid  at  300  MHz: 

Sll  =  S22  =  -0.1191  +  j  0.0637 

SI  2  =  S21  =  0 

EAl  1  =  EA22  =  0.229 

EAl  2  =  EA2 1  =  0 

Although  the  piecewise-linear  integration  technique  is  efficient 
and  successful  for  most  problems,  it  is  not  appropriate  in  forward 
scattering  problems  (where  w  =  0  everywhere  on  the  target)  or  in 
specular  scattering  from  a  flat  region  on  the  target  (where  w  is  constant 
over  a  significant  portion  of  the  surface).  Therefore,  the  computer  out¬ 
put  data  are  not  reliable  in  these  cases. 


VII  .  COMPUTING  THE  PULSE  RESPONSE 

The  pulse  scattering  problem  is  treated  in  the  fourth  and  last 
section  of  the  program  shown  in  Fig.  7  .  The  entire  program  is  ready 
for  the  computer  when  the  four  sections  are  put  together  (with  the  card 
numbers  running  sequentially  from  1  through  345)  and  backed  up  with 
suitable  input  data  (such  as  cards  1  through  299  in  Fig.  1). 

To  explain  the  pulse  calculations,  it  is  convenient  to  write  Eq.  (49) 
in  the  following  form 
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3  READ (5.29 >FMC »L»M  212 

WAVE=299.79/FMC  213 

WR  I  TE  (  6 . 32  )  FMC  »  WAVE  214 

TPL=TP/WAVE  215 

SX  = ( .O. .0 )  216 

SY=(.C..O)  217 

SZ=(.0,.0>  218 

GK=TPL*W1  219 

GA  =  TPL*  WM 1 N  220 

DG=  TPL  *DEL  22 1 

0 A  =  CMPLX ( COS ( GA ) . S I N ( GA )  )  222 

DO  240  K  = 1 »KX  223 

G8=GA+DG  224 

QB=CMPLX (COS ( G8 ) . S I N (GB ) )  225 

GL  =GK+DG  226 

OC*GB*OB-GA*OA  227 

OD=OB-OA  228 

OK=QU*GL-GC-( .0. I • )*QD  229 

QL*OC-OD*GK+ ( .0. 1 . ) *OD  230 

SX*  SX  +  QK*  AX  ( K  ) +QL*  AX  (  K-f  1 )  231 

SY*SY+QK*AY(K)+QL*AY(K+1 )  232 

SZ=SZ+QK*AZ (K )+QL*AZ(K+l )  233 

GA  =G8  234 

GK*GL  235 

240  OA  =  08  236 

DEN=TP*DG  237 

S 1 1  * (SX*X1 1 +SY*Y1 1 ♦SZ*Zl 1 ) /DEN  238 

SI  2= (SX*X12  +  SY*Y12  +  SZ*Z12 ) /DEN  239 

S21 =SZ*Z2l/DEN  240 

S 22= (SX*X22+SY#Y22+SZ*Z22 )/DEN  241 

WR 1 TE ( 6 *  33 )  SI  1 . S 1 2 . S2 1  * S22  242 

EA1 1 = (CABS(S1 1 >**2 >*FPl  243 

E A  1 2= ( CABS (S12)**2)*FP1  244 

EA2 1 = ( CABS ( S2 1 )**2)*FPl  245 

EA22* (CABS(S22 )**2 )*FP1  246 

WR 1 TE ( 6 . 34 )  EA1 1 .EA12.EA21 .EA22  247 

WRITE  (6.23)  248 


Fig.  6.  Third  section  of  the  computer  program. 
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(68) 


F(t)  =  _L 

2tt 

where 

(69)  g  =  kw  , 

(70)  xa  =  GN  or  -cot  (whichever  is  larger)  , 

(7  1)  y a  =  GX  or  ojt  -  cot  (whichever  is  smaller)  , 

(7  2)  GN  =  k  WMIN, 

and 

(73)  GX  =  k  WMAX  . 


pVa 

rYa 

cos  ojt  \  A(g)  cos  g  dg  -  sinwt 

— 

\  A(g)  sing  dg 

xa 

The  area  function  A(g)  is  the  same  as  A(w),  and  it  has  already  been 
calculated  and  stored  in  AX(K),  AY(K)  and  AZ(K). 

The  integer  L,  determines  the  number  of  points  per  cycle  to  be 
calculated  on  the  scattering  waveforms.  Card  254  calculates  the  size 
of  the  increments  in  cot  as  follows: 


(74)  DWT  =  2tt/l  . 


The  integer  M  specifies  the  number  of  complete  cycles  in  the  incident 
waveform.  The  pulse  width  of  the  incident  waveform,  as  measured  in 
radians  and  observed  at  an  arbitrary  point  in  space,  is  calculated  by 
card  253  as  follows: 


(75)  coT  =  WTAU  =  2ttM 
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The  scattering  function  F(t)  is  nonzero  only  in  the  following  interval: 
-  GX  <  cot  <  co  t  -  GN  • 

It  is  convenient  to  assign  symbols  to  the  integrals  in  Eq.  (68) 
as  follows: 


and 


pY  a 

(77)  V  (t)  =  \  A(g)  sing  dg  . 

xa 

These  are  functions  of  time  by  virtue  of  the  time  dependent  limits  xa 
and  ya*  In  the  program  the  rectangular  components  of  U(t)  are  denoted 
by  UX,  UY  and  UZ,  and  the  components  of  V(t)  are  VX,  VY  and  VZ. 

WT,  SWT  and  CWT  denote  cot,  sincot  and  coscot. 

U(t)  and  V(t)  are  identically  zero  for  all  values  of  cot  less  than 
or  equal  to  -  GX .  Thus,  the  computer  begins  the  pulse-response 
calculations  by  setting  cot  =  -  GX  (card  261),  initializing  U(t)  and  V(t) 
to  zero  (cards  269  -  274)  and  writing  out  the  first  data  point  on  the  scat¬ 
tering  waveform  (card  268).  At  this  point  in  time,  the  integration  limits 
are  found  from  Eqs  .  (7  0)  and  (7  1)  to  be  xa  =  ya  =  GX . 

Suppose  the  computer  has  just  completed  the  calculation  of  the 
scattering  data  for  a  given  value  of  cot,  has  incremented  cot  in  card  277, 
and  is  preparing  to  calculate  the  next  point  on  the  scattering  waveform. 
The  first  step  in  the  procedure  is  to  determine  the  integration  limits 
xa  and  ya,  denoted  by  XA  and  YA.  This  is  accomplished  in  cards  281  - 
292  in  accordance  with  Eqs.  (70)  and  (71).  These  cards  also  determine 
the  indices  KA  and  LA  of  the  segments  on  the  g  axis  which  contain  the 
points  g  =  XA  and  g  =  YA,  respectively. 

The  integrals  U(t')  and  V(t!)  currently  in  storage  have  limits  now 
denoted  by  XB  and  YB  by  virtue  of  cards  338  and  339*  Since  XA  and  YA 
are  equal  to  or  less  than  XB  and  YB,  respectively, 


(78)  U(t)  =  U(t')  + 
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and 


pXB  pYB 

(79)  V(t)  =  V(t,)  +  \  A(g)  sin  g  dg  -  \  A(g)singdg 

JXA  JYA 


Calculations  based  on  Eqs  .  (78)  and  (79)  are,  of  course,  far  more  efficient 
than  Eqs.  (76)  and  (77). 

Equations  (78)  and  (79)  are  programmed  in  cards  294  -  329-  The 
piecewise -linear  technique  is  employed  to  evaluate  the  integrals  in  these 
equations.  The  integrals  with  limits  XA  and  XB  are  evaluated  in  the 
first  pass  through  the  DO  LOOP  that  begins  with  card  299,  and  those 
with  limits  YA  and  YB  are  calculated  in  the  second  pass.  In  each  pass, 

ZA  and  ZB  denote  the  limits  of  integration.  Initially  I  denotes  the  index 
number  of  the  g-axis  segment  containing  the  point  g  =  ZA  and  G  is  the 
value  of  g  at  the  center  of  this  segment.  A(g)  is  treated  as  a  linear 
function  over  this  segment  in  the  same  manner  as  in  Eqs.  (55)  through 
(57).  GP  denotes  the  value  of  g  at  the  center  of  segment  I  +  1,  and  HA 
and  HB  are  the  limits  of  the  subintegral  on  segment  I.  When  the  inte¬ 
gration  on  segment  I  is  finished,  the  computer  moves  on  to  the  next 
segment  (unless  ZA  and  ZB  lie  on  the  same  segment).  The  following 
equations  are  used  in  cards  310  -  319  to  evaluate  the  subintegral  on 
segment  I: 

(80)  (  ^(g)  cos  g  dg  =  _L  [  Cl  A(I)  +  CII  A(I+1)]  , 

~ha~  dg 

(81)  \  A(g)  sin  g  dg  =  — L_  [  DI  A(I)  +  DII  A(I+1)J  , 

Jh  DG 

na 


(82)  Cl  =  (ha  -  gi+i)  sinha  +  (gi+i  -  hb)  sinhb  +  cos  h  a  -  cos  hb  , 

(83)  CII  =  (gi  -  ha)  sinha  +  (hb  -  g^)  sinhb  +  cos  hb  -  cos  ha  , 


(84)  DI  =  (gi+i-  ha)  cos  ha  +  (hb  -  gi+i)  cos  hb  +  sinha  -  sinhb  , 


and 
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(85) 


DII  =  (ha  -  gi)  cosha  +  (gi  -  hb)  cos  hb  +  sinhb  -  sinha  . 


DG  is  defined  in  Eq .  (59)-  In  the  program,  g{,  gi+i  ,  ha  and  hb  are 
denoted  by  G,  GP,  HA  and  HB  .  SGN  is  positive  in  the  first  pass  through 
this  integration  loop  and  negative  in  the  second  pass,  in  accordance  with 
the  plus  and  minus  signs  in  Eqs  .  (78)  and  (79)  • 

When  the  vectors  U(t)  and  V(t)  have  been  calculated,  Eq .  (68)  is 
used  to  determine  F(t)  .  The  rectangular  components  of  F(t)  are  denoted 
by  FX,  FY  and  FZ  in  cards  330-332.  Finally  cards  333  -  336  calculate 
the  scattering  functions  defined  by  Eqs.  (43)  through  (47),  and  card  337 
arranges  the  writeout  of  these  data  on  the  scattering  waveforms  for  the 
pulse  situation. 

The  results  obtained  in  this  manner  (using  the  input  data  in  Fig.  1) 
for  a  prolate  spheroid  are  illustrated  in  Fig.  8. 


VIII.  CONCLUSIONS 

A  digital  computer  program  is  described  which  uses  the  physical 
optics  approximation  for  scattering  by  convex  perfectly  conducting  tar¬ 
gets  with  arbitrary  shape.  The  target  size  and  shape  are  specified  by 
input  data  giving  the  coordinates  of  many  points  on  the  surface.  This 
program  handles  bistatic  and  backscatter  problems.  The  input  data 
specify  the  frequency,  the  incidence  angles  ( 0 ,  4>i )  and  the  scattering 
angles  (0s*  ^s)*  In  the  CW  case,  the  output  data  give  the  complex 
elements  in  the  scattering  matrix.  Typical  results  are  shown  in  the 
Appendix  for  scattering  from  spheres,  spheroids,  ogives  and  cones. 

This  scattering  program  also  handles  the  pulse  case  where  the 
incident  wave  has  a  finite  number  of  cycles.  For  this  case,  the 
following  additional  input  data  are  required:  the  number  of  cycles  in 
the  incident  wavetrain,  and  the  number  of  points  per  cycle  to  be  computed 
for  the  scattering  waveforms.  The  output  data  give  the  four  waveforms 
associated  with  the  four  elements  in  the  scattering  matrix. 

Simple  modifications  can  be  made  in  the  program,  if  desired,  to 
increment  the  frequency,  the  incidence  angles  or  the  scattering  angles  . 
This  is  convenient  for  generating  various  types  of  scattering  curves 
such  as  those  in  the  Appendix. 
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4  CONTINUE  24 9 

LM=L*M  250 

FL=L  251 

FM=M  252 

WTAU=TP*FM  253 

DWT=TP/FL  254 

IF (LM.tU.« >G0  TO  300  255 

DGT  =  TPL*DELT  256 

G1=TPL*WI  257 

GX=TPL*WMAX  258 

GN=TPL*WM1N  259 

DEN=TP#DG  260 

WT  A  =  — GX  261 

WTB=WTAU-GN  262 

*T=WTA  263 

FI  1=0.0  264 

F 1 2=  0 • 0  265 

F2 1 =0.0  266 

F  22  =  0 • 0  267 

WR1TE<6.22)  'JlT.Fl  1  .F12.F21  .F22  268 

UX  =  .  0  269 

UY=.0  270 

UZ=.0  271 

VX= . 0  272 

VY=.0  273 

VZ=.0  274 

X9=GX  275 

YB=GX  276 

260  nt/T  =  WT  +OWT  277 

IF  <*/T  *GT  •  wTt3  )  GO  TO  300  278 

SWT=SIN<WT>  279 

CWT  =COS  <  WT  >  280 

X A=GN  281 

KA= 1  282 

XX  =  —  WT  283 

IF (XX.LE.XA)  GO  TO  265  284 

X  A  =  XX  285 

KA= (XA-G1  J/DG+l .5  286 

265  Y A  =  GX  287 

LA=KX  288 

YY=WTAU-WT  289 

lF(YY.GE.YA)  GO  TO  268  290 

Y  A  =  Y Y  291 

l_A=  <  YA-G  1  l/DG+l  .5  292 

268  SGN=1.  293 

Z A  =  XA  294 

ZB  =  XB  295 

1 =KA  296 

F  =  I  —  1  297 

G  =  F  *DG+-G  1  298 

DO  290  K=l .2  299 

lF(ZA.GE.ZB)  GO  TO  285  300 

270  GP=G+DG  301 

HA  =  G~DGT  302 

HB=G+DGT  303 

IF(HA.LT.ZA)  HA=ZA  304 

1  F  (Hd.GT.Zo  )  H6  =  Zt)  305 

SA=S1N(HA)  306 

SB  =  S 1 N ( HB  >  307 

C A  =COS ( HA  >  308 

CB  =  COS ( HB  >  J09 

C  1  =  <  SA*  (HA-GP  >  +  St3*  *  GP-HB  )  +  C  A  —  CB  >  *SGN  310 

Fig.  7  .  Last  section  of  the  computer  program. 
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C  I  I  =  (  SA*  (G-HA  >+SU*  (Md-Ol+Cu-CA  >  «SGN  Jt  1 

UX=UX+CI*AX<  1  )+CI  1*AX(  1  +  1  )  312 

UY=UY+C I *AY (  I  )  +C I  I *AV (  I  +  1  )  313 

UZ=UZ+C  I  #AZ  (  I  )+ClI*AZ(I+|  )  314 

D  I  =  (  SA-SB  +  CA*  (GP-HA  >+CB*  (Hb-GP  )  )  *SGN  31  5 

D  1  1  -  (  SB— SA+CA*  (MA-G  )  +  CB*  (G-HiJ  )  )  *SGN  316 

VX  =  VX  +  D1*AX(  I  ) +D I  I  * AX (  1  +  1  )  317 

VY  =  VY  +  D1*AY<  1  >+DI  I  *AY  <  1  +  1  )  318 

VZ=VZ+DI *AZ ( I ) +DI 1 *AZ ( I + I )  319 

lFlHd.GE. ZB)  GO  TO  28b  320 

1*1+1  321 

G=GP  322 

GO  TO  270  323 

285  ZA=YA  324 

ZB=YB  325 

1=LA  326 

F  = 1 -  1  32  7 

G  =  F  *DG  +  G I  328 

290  SGN=- 1 .  329 

FX= <UX*CWT-VX*SWT )/DEN  330 

FY= <UY*CWT-VY*SWT J/DEN  331 

FZ= <UZ*CWT-VZ*SWT J/DEN  332 

F 1 1 =  FX*X 1 1 +FY*Y 1 1 +FZ*Z 1 1  333 

F 1 2*FX*X 1 2+FY*Y 1 2  +  FZ*Z  1  2  334 

F  2 1 =FZ*Z2 1  335 

F22=FX*X22+FY#Y22+FZ*Z22  336 

WR 1 T  E ( 6 «  22  >  WT.F11  .F12.F2I  .F22  337 

XB  =  X  A  338 

YB=YA  339 

GO  TO  260  340 

300  CONTINUE  341 

WRITE  (6.23)  342 

STOP  343 

END  344 

SDAT  A  345 


Fig.  7.  Last  section  of  the  computer  program,  (cont.) 
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Fig.  8.  Calculated  pulse  response  for  prolate  spheroid, 
using  input  data  from  Fig.  1. 
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APPENDIX 

NUMERICAL  RESULTS 


This  Appendix  presents  the  results  obtained  with  the  physical- 
optics  computer  program  for  several  target  shapes.  Table  I  lists  the 
number  of  planes  used  to  describe  each  target,  the  total  number  of 
points  on  the  target,  and  the  time  required  for  an  IBM  7094  to  process 
the  target  description  data  and  set  up  the  polyhedral  approximation. 
The  last  column  lists  the  additional  computer  time  needed  to  calculate 
the  area  function  and  the  CW  scattering  matrix. 

TABLE  I 

STATISTICS  FOR  THE  TEST  CASES 


Target  Shape 

Planes 

Points 

Time  for 
Polyhedron 

Area  Function 

and  CW  Matrix 

Sphere 

30 

1083 

8 . 6  sec 

1  . 6  sec 

Prolate 

Spheroid 

30 

774 

6.4 

1.3 

Oblate 

Spheriod 

30 

1257 

9.9 

1  -9 

Ogive 

50 

1008 

8.2 

1  .8 

Cone 

35 

1176 

9.8 

3.3 

The  dimensions  of  the  targets  are  specified  in  Figs.  9  through  13. 

It  is  found  from  Table  I  that  the  computation  time  for  processing 
the  target  description  data  and  setting  up  the  polyhedron  increases 
linearly  with  the  number  of  points  on  the  target  surface.  The  computer 
handles  about  125  points  (to  establish  250  facets)  per  second,  regardless 
of  the  shape  of  the  target. 

In  generating  the  area  function  and  the  CW  scattering  matrix,  the 
computer  spends  very  little  time  on  the  shadowed  facets.  Thus,  the 
computation  time  for  this  work  increases  linearly  with  the  number  of 
illuminated  facets  on  the  target.  This  explains  why  the  computer  spent 
twice  as  much  time  on  the  cone  (3.3  seconds)  as  on  the  sphere  (1  .6  seconds) 
as  shown  in  Table  I  This  sphere  has  approximately  2000  facets,  and  half 
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of  them  are  illuminated.  The  cone  also  has  approximately  2000  facets, 
and  all  of  them  are  illuminated  for  axial  incidence.* 

In  calculating  the  area  function,  the  computer  handles  about  670 
illuminated  facets  per  second. 

To  use  the  computer  program  most  efficiently,  one  would  cal¬ 
culate  a  significant  amount  of  scattering  data  for  one  target  before 
passing  on  to  the  next  one.  Furthermore,  to  minimize  the  are  a- function 
calculations,  one  would  calculate  the  scattering  data  for  many  different 
frequencies  before  changing  0^,  cj>{,  0S  or  c(>s  •  When  these  incidence 
and  scattering  angles  are  held  constant,  the  computer  generates  CW 
scattering  data  for  20  different  frequencies  in  one  second  for  the  oblate 
spheroid  in  Table  I,  and  for  10  frequencies  per  second  for  the  cone. 

For  the  sphere,  it  generates  CW  and  pulse  -  response  data  for  a  new 
frequency  in  one  second  (with  three  cycles  in  the  incident  wavetrain 
and  20  calculated  points  per  cycle  in  the  scattering  waveform). 

Figures  9  through  13  display  the  numerical  results  for  scattering 
by  several  different  target  shapes.  The  dots  represent  the  output  data 
from  the  computer  program  described  herein.  The  input  data  listed 
the  coordinates  of  a  large  number  of  points  (Table  I)  distributed  almost 
uniformly  over  the  target  surface.  The  source  of  the  physical-optics 
data  (the  solid  curves  in  Figs.  9  -  13)  is  described  a  little  farther  on 
in  this  Appendix. 

Figures  9a  through  13a  show  curves  of  backscatter  echo  area 
(in  square  meters)  versus  frequency  for  axial  incidence.  For  each  of 
these  targets,  it  may  be  noted  that  the  computer  output  agrees  closely 
with  the  physical  optics  solution  for  all  frequencies  up  to  600  MHz. 

(This  simply  indicates  that  the  computer  program  does  a  good  job  in 
generating  physical  optics  data.  It  is  not  implied  that  the  physical 
optics  formulation  always  gives  a  good  approximation  to  the  exact 
scattering  data.)  To  obtain  equally  good  results  at  higher  frequencies, 
it  would  be  necessary  to  include  more  points  in  the  input  data  description 
of  the  target . 


*  For  all  the  targets  listed  in  Table  I  except  the  cone,  the  input  data 
points  covered  the  whole  surface  and  the  computer  set  up  a  closed  poly¬ 
hedron.  For  the  cone,  however,  the  input  data  did  not  cover  the  flat 
plate  at  the  base.  Since  the  base  is  shadowed  for  all  cases  shown  in 
Fig.  13,  the  computed  data  apply  to  a  cone  with  a  circular  disk  at  the  base. 
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Although  the  computer  program  has  not  been  adequately  tested 
for  the  most  general  target  shapes  and  bistatic  situations  for  which  it 
is  designed,  it  is  believed  that  it  will  give  accurate  physical-optics 
data  except  for  concave  targets,  forward  scattering,  and  specular 
scattering  from  a  large  flat  subsurface.  The  program  could  be 
modified  without  much  difficulty  to  remove  the  last  two  of  these 
exceptions  . 

All  of  the  targets  considered  in  this  Appendix  have  axial  symmetry 
with  respect  to  the  z  axis.  To  check  the  accuracy  of  the  computer  out¬ 
put  data  shown  in  Figs.  9a  -  13a,  closed-form  expressions  were  derived 
for  the  vector  area  function  for  each  of  these  targets  for  the  backscatter 
situations  with  axial  incidence.  Letting  =  4>i  =  6S  “  <j>s  =  0,  we  find 
from  Eq.  (22)  that  w  =  2z  .  From  symmetry  considerations,  Ax(w )  - 
Ay(w)  =  0.  For  targets  having  axial  symmetry  with  respect  to  the  z 
axis,  it  is  convenient  to  describe  the  target  shape  with  an  equation  of 
the  form  p  =  p(z)  where  (p,cj>,  z)  represent  the  coordinates  in  the 
cylindrical  system.  Then  it  is  easy  to  show,  with  the  aid  of  Eq.  (51), 
that 

(84)  Az(w>  =  -  Trp  p  '  , 


where  p1  represents  the  derivative  of  p(z)  with  respect  to  z.  Equation  (84) 
applies  over  the  illuminated  region  of  the  target  where  p1  is  negative, 
and  Az(w)  =  0  over  the  shadowed  region  where  p'  is  positive. 

All  the  targets  considered  in  this  Appendix  can  be  represented 
with  the  following  equation: 


(85)  F(p,  z)  =  A^2  +  A2z2  +  A3  pz  +  A*  z  +  A5p  +  A6  =0  . 


From  Eqs  .  (84)  and  (85),  the  area  function  for  axial  backscatter  is 
given  in  closed  form  by 


(86)  Az  =  tt  p  ( 2 A2 z  +  A3p  +  A4)/(2Aip  +  A3z  +  A5  )  . 


For  spheres,  prolate  spheroids  and  oblate  spheroids, 
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(87) 


Az  ~ 


TTZ  B2 

~7F~ 


where  A  and  B  denote  the  semimajor  and  semiminor  axes,  respectively. 
For  the  ogive, 

(88)  Az  =  T.PZ„. 

vjB2  -  z2 

where  B  is  the  radius  of  the  circular  are  that  generates  the  ogive. 

For  the  cone, 


(89)  Az  =  tt  z  tan2  cx  , 


where  0C-  denotes  the  half- angle  of  the  cone. 

For  axial  backscatter  from  a  target  with  axial  symmetry,  it  is 
found  from  Eq .  (54)  that  Sx  =  Sy  =  0  and 

(90)  Sz  =  (2j/\)  (  Az  eZjkz  dz  . 


The  integral  in  Eq .  (90)  can  be  evaluated  to  obtain  closed-form  expressions 
for  the  scattering  matrix  and  the  echo  area  for  the  cone  and  the  spheroids. 
The  resulting  data  are  shown  by  the  solid  curves  with  the  label  "physical 
optics”  in  Figs.  9a,  10a,  11a  and  13a.  For  the  ogive  in  Fig.  12a,  the 
"physical  optics"  curve  was  obtained  from  Eqs  .  (88)  and  (90)  with 
numerical  integration. 

Figure  9b  shows  backscatter  results  for  the  sphere  as  a  function 
of  the  aspect  angle.  For  a  true  sphere,  the  echo  area  is  of  course 
independent  of  aspect  angle.  Thus,  the  variations  in  the  computer  out¬ 
put  data  arise  from  the  polyhedral  approximation. 

Figures  1  Ob  -  13b  show  the  backscatter  results  as  a  function  of  the 
aspect  angle  (0^  =  0g  =  0)  for  the  spheroids,  ogive  and  cone.  The  solid 
curves  with  the  label  "physical  optics11  were  obtained  with  another  com¬ 
puter  program  which  is  designed  especially  for  backscattering  from  tar¬ 
gets  with  axial  symmetry.  The  input  data  for  this  auxiliary  program 
describes  the  target  shape  by  assigning  numerical  values  to  the  coefficients 
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in  Eq.  (85).  The  angles  c^i  and  <\>s  are  assumed  to  be  zero,  the  cylindri¬ 
cal  coordinates  (p,  <j),  z)  are  employed  for  points  on  the  target  surface, 
and 


(91)  w  =  2(p  sin0  cos(j)  +  z  cos0). 


It  is  found  that  Sy  =  0, 


(92)  S  =  (jA)  j  j  (x  COS  4)  -  z  P ' )  P  eJkw  d(j)  dz  , 


(93)  Sn  =  S22  =  -  Sx  sin 0  -  Sz  cos  0  , 


and  S 12  =  S21  =  0.  The  auxiliary  program  uses  numerical  integration 
to  evaluate  S  in  Eq.  (92). 

Bistatic  scattering  results  are  shown  in  Figs.  9c  -  13c  and  9d  - 
13d.  These  curves  apply  for  axial  incidence  with  0i  -  <j>i  =  ^s  =  so 
that 

(94)  w  =  p  sin0s  cos  cj)  +  z(  1  +  cos  0S)  , 

f  ikz(  1  +  cos  0  o) 

(95)  Sx  =  -  k  \  p  J i(k  p  sin  0  s )  e  dz  , 

(96)  Sy  =  S12  =  S21  =  o  , 

C  .  ikz(l  +  cos0c) 

(97)  Sz  =  -  jk  J  p  p*  J Q(k  p  sin  0S )  eJ  dz  , 

(98)  Sn  =  -  Sx  sin  0  s  -  Sz  cos  9S  , 

(99)  S22  =  -  Sz  , 

and  JQ  and  J*  represent  Bessel  functions. 
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Equations  (94)  -  (99)  apply  only  for  axial  incidence  on  a  target  having 
axial  symmetry.  (Similar  expressions  are  given  in  Reference  2.)  These 
equations  were  used  to  develop  another  auxiliary  computer  program  to 
generate  the  E-plane  and  H-plane  curves  with  the  label  ’’physical  optics”. 
This  program  also  used  the  coefficients  in  Eq.  (85)  as  input  data  for  the 
target  shape,  and  evaluated  Sx  and  Sz  with  numerical  integration. 

The  bistatic  curves  in  Figs.  9-13  agree  with  the  physical -optics 
echo  area  expression  for  the  forward- scattering  situation: 

(100)  0-  =  4ttA2/\2  , 


where  A  represents  the  area  of  the  target  projected  along  the  line  of 
sight  between  the  transmitting  and  receiving  antennas.  In  the  forward¬ 
scattering  situation  (and  also  for  backsc atte ring) ,  the  physical  optics 
echo  area  shows  no  polarization  dependence  even  for  a  target  with 
arbitrary  shape.  Equation  (100)  is  not  restricted  to  symmetrical  tar¬ 
gets  or  axial  incidence. 

Though  all  the  data  in  Figs.  9-13  were  obtained  with  the  aid  of 
various  digital  computer  programs  based  on  the  physical-optics  con¬ 
cepts,  the  close  agreement  between  the  dots  (with  the  label  "computer 
output”)  and  the  curves  (with  the  label  ’’physical  optics”)  is  significant 
for  two  reasons  .  First,  the  general  program  uses  a  target  description 
technique  which  differs  greatly  from  that  used  in  the  specialized  programs  . 
Furthermore,  some  of  the  specialized  programs  used  closed-form  ex¬ 
pressions  and  others  were  based  on  equations  involving  line  integrals  or 
surface  integrals,  whereas  the  general  program  handles  each  problem  in 
the  same  manner  (setting  up  the  polyhedron,  calculating  the  vector  area 
function,  and  integrating  the  area  function  to  obtain  the  scattering  mat¬ 
rix)  .  Having  demonstrated  that  this  general  procedure  gives  correct 
physical- optics  data  for  a  variety  of  target  shapes  and  bistatic  and 
backscatter  situations,  we  have  confidence  that  the  program  will  perform 
equally  well  with  other  targets  and  situations. 
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ECHO  AREA  a..  -  ECHO  AREA 


(a)  BACKSCATTER  vs.  FREQUENCY  (b)  BACKSCATTER  vs.  ASPECT  AT  300  MHz 


(c)  E-PLANE  BISTATIC  SCATTERING  AT  300MHz 


(d)  h-PLANE  BISTATIC  SCATTERING 
AT  300  MHZ 


Fig.  9.  Scattering  data  for  a  sphere  with  a  radius  of  one  meter. 
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ECHO  AREA  <r,,  ECHO  AREA 


(a)  BACKSCATTER  vs.  FREQUENCY 
FOR  AXIAL  INCIDENCE 


(b)  BACKSCATTER  vs.  8  AT  300  MHz 


(C)  E-  PLANE  BISTATIC  SCATTERING 
AT  300  MHz  FOR  AXIAL  INCIDENCE 


(d)  H- PLANE  BISTATIC  SCATTERING 
AT  300  MHz  FOR  AXIAL  INCIDENCE 


Fig.  10. 


Scattering  data  for  a  prolate  spheroid 
with  A  =  1  m  and  B  =  0.5  m. 


39 


ECHO  AREA  <r..  -  ECHO  AREA 


(a)  8ACKSCATTER  vs  FREQUENCY  FOR  AXIAL  INCIDENCE  < b)  BACKSCATTER  vs  9  AT  300  MHz 


0  60  120  100  240 

8%  (DEGREES) 

(C)  E-PLANE  BISTATIC  SCATTERING  AT  300  MHz 
FOR  AXIAL  INCIDENCE 


0,  (DEGREES) 

(d)  h- plane  bistatic  scattering  at 

300  MHz  FOR  AXIAL  INCIDENCE 


Fig.  11  .  Scattering  data  for  an  oblate  spheroid 
with  A  =  0.5m  ana  B  =  1  m . 
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ECHO  AREA  <ru  P  ECHO  AREA 


(a)  BACKSCATTER  vt  FREQUENCY  FOR  AXIAL  INCIDENCE 


(C)  E-  PLANE  BlSTATtC  SCATTER  NG  AT  300  Mhi 
FOR  AXIAL  INCIOENCE 


(d)  H  — PLANE  B1STAT1C  SCATTERING  AT  300MHz 
FOR  AXIAL  INCIDENCE 


Fig.  12.  Scattering  data  for  an  ogive  with 
A  =  1  m  and  ol  =  22.5°. 
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ECHO  AREA 


(a)  BACKSCATTER  vs  FREQUENCY 
FOR  AXIAL  INCIDENCE 


(c)  E-  PLANE  BISTATIC  SCATTERING  AT 
300  MHz  FOR  AXIAL  INCIDENCE 


0  30  60  90 

9  (  DEGREES  ) 


(b)  BACKSCATTER  vs  9  AT  300  MHz 


( d  )  H- PLANE  BISTATIC  SCATTERING  AT 
300  MHz  FOR  AXIAL  INCIDENCE 


Fig.  13.  Scattering  data  for  a  cone  with 
L  =  1  m  and  cx  =  18°  . 
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